06 - Multiple lineare Regression

Auswertung Emprischer Daten

Prof. Dr. David F. Urschler

Vorgehensweise

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

Welche Fragen können wir beantworten?

  • Welche Faktoren sagen am besten den Resozialisierungserfolg vorher?

  • Welche Faktoren haben einen Einfluss auf Prdoduktbewertungen?

  • Welche “Umstände” beinflussen einen Krankheitsverlauf?

  • Was sind die besten Prädiktoren für Studienerfolg?

  • Was sind Indikatoren für Unternehmenserfolg?

Beispiel

Werden Menschen aggressiver, wenn …

  • sie Helene Fischer häufiger hören?
  • sie eine Abneigung verspüren?

Information zum Datensatz

ae_mult_reg.csv

  • demographische Daten (id, age & gen: 0 = weiblich, 1 = männlich)
  • Häufigkeit der Exposition an einem Abend (hexp: 0 bis 15 Mal an einem Abend)
  • Abneigung gegen den Song (relu: 0 = keine Ablehnung, 3 = ich springe gleich aus dem Fenster)
  • Menge an scharfer Sauce in Gramm (mghs: 0 bis 35 Gramm)

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style A fill: #009e73

Daten einlesen

Zeig mir den Code
ds <- read.csv(here::here("ae_mult_reg.csv"),
               header = TRUE)

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style B fill: #009e73

Übersicht verschaffen



id age gen hexp relu mghs
1 22 1 15 2.1 26.91
2 40 0 8 1.8 17.85
3 39 0 10 1.9 29.37
4 42 1 10 1.1 12.01
5 38 1 7 0.7 14.31
6 21 1 3 1.1 9.67

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style C fill: #009e73

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style D fill: #009e73

Ergeben die Daten Sinn?

Prädiktoren

Zeig mir den Code
# Häufigkeit
summary(ds$hexp)

# Abneigung
summary(ds$relu)

Ergeben die Daten Sinn?

Kriterium

Zeig mir den Code
# Menge an scharfer Sauce
summary(ds$mghs)

Ergeben die Daten Sinn?

demographische Daten

Zeig mir den Code
# Gender
table(ds$gen)

# Alter
summary(ds$age)

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style E fill: #009e73

Skalenwerte berechnen

Demographische Variablen als Faktor

Zeig mir den Code
# Gender
ds$gen_f <- factor(ds$gen,
                   levels = c(0, 1),
                   labels = c("w", "m"))

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style F fill: #009e73

Multiple lineare Regression berechnen

Voraussetzungen

  • Modell ist richtig spezifiziert
  • Lineare Beziehung zwischen den Variablen
  • Normalverteilung der Residuen
  • Homoskedastizität der Residuen
  • Unabhängigkeit der Residuen (Autokorrelation)
  • keine Multikollinearität
  • keine Ausreisßer

Multiple lineare Regression berechnen

Voraussetzungen

Wir berechnen zuerst das Modell und prüfen die Vorrausetzungen im Anschluss

Multiple lineare Regression berechnen

Berechnen des Modells

Zeig mir den Code
m1 <- lm(mghs ~ hexp + relu, # "+" ohne Interaktion
         data = ds)

Multiple lineare Regression berechnen

Berechnen des Modells

Was haben wir gefunden?

Zeig mir den Code
summary(m1)

Multiple lineare Regression berechnen

Berechnen des Modells

Was haben wir gefunden


Call:
lm(formula = mghs ~ hexp + relu, data = ds)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.239  -3.413  -0.077   3.705  10.560 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.42357    0.77187  16.095   <2e-16 ***
hexp        -0.02665    0.08974  -0.297    0.767    
relu         5.40432    0.57024   9.477   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.655 on 272 degrees of freedom
Multiple R-squared:  0.2881,    Adjusted R-squared:  0.2829 
F-statistic: 55.04 on 2 and 272 DF,  p-value: < 2.2e-16

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Graphisch

Zeig mir den Code
plot(m1, 1)

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Statistisch I - RESET

  • steht für: Ramsey Regression Equation Specification Error Test
  • fügt Prädiktoren höhrere Ordnung hinzu
  • resettest()-Funktion aus dem lmtest-Package

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Statistisch I - RESET

Zeig mir den Code
lmtest::resettest(m1,
                  power = 2:3,  # fügt quadratische und kubische Terme hinzu
                  type = "regressor") # prüft für die Prädiktoren

    RESET test

data:  m1
RESET = 12.84, df1 = 4, df2 = 268, p-value = 1.411e-09

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Statistisch 2 - Rainbow-Test

  • prüft für mittlere Ausprägungen der Prädiktoren
  • raintest()-Funktion aus dem lmtest-Package

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Statistisch 2 - Rainbow-Test

Zeig mir den Code
lmtest::raintest(m1)

    Rainbow test

data:  m1
Rain = 1.2102, df1 = 138, df2 = 134, p-value = 0.1339

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Vorrausetzun scheint verletzt zu sein, was nun?

  • Anderen Voraussetzungen prüfen
  • Anpassen des Modells durch hinzüfügen/entfernen von Prädiktoren bzw Interaktionen
  • Polynominale Terme

  • … und nochmal

Voraussetzungen

Normalverteilung der Residuen

Graphisch - Histogramm der Residuen

Zeig mir den Code
hist(m1$residuals)

Voraussetzungen

Normalverteilung der Residuen

Graphisch - Q-Q Plot

Zeig mir den Code
plot(m1, 2)

Voraussetzungen

Normalverteilung der Residuen

Statistisch - Shapiro-Wilk-Test

Zeig mir den Code
shapiro.test(m1$residuals)

    Shapiro-Wilk normality test

data:  m1$residuals
W = 0.98724, p-value = 0.01545

Voraussetzungen

Homoskedastizität

Graphisch

Zeig mir den Code
plot(m1, 3)

Voraussetzungen

Homoskedastizität

Statistisch - Breusch-Pagan-Test

  • ncvTest()-Funktion aus dem car-Package

Voraussetzungen

Homoskedastizität

Statistisch - Breusch-Pagan-Test

Zeig mir den Code
car::ncvTest(m1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.29104, Df = 1, p = 0.25586

Voraussetzungen

Unabhängigkeit der Residuen (Autokorrelation)

  • nur bei Längsschnittdesigns und gelcusterten Daten von Relevant
  • Durbin-Watson-Test

Voraussetzungen

Statistisch

Zeig mir den Code
lmtest::dwtest(m1)

    Durbin-Watson test

data:  m1
DW = 1.9103, p-value = 0.2279
alternative hypothesis: true autocorrelation is greater than 0

Voraussetzungen

Multikollinearität

  • Idee, dass jeder Prädiktor einen eigenständigen Beitrag hat

Voraussetzungen

Multikollinearität

Korrelationsmatrix

Variable M SD X1 X2
1. hexp 7.59 3.52
2. relu 1.22 0.55 .45**
3. mghs 18.81 5.50 .23** .54**

Voraussetzungen

Multikollinearität

Variance inflation factor (VIF)

  • ncvTest()-Funktion aus dem car-Package
  • VIF ≈ 1 -> Prädiktoren sind nicht korreliert
  • VIF ≈ 5 -> Prädiktoren sind moderat korreliert
  • VIF ≈ 10 -> Prädiktoren sind hoch korreliert

Voraussetzungen

Multikollinearität

Variance inflation factor (VIF)

Zeig mir den Code
car::vif(m1)
   hexp    relu 
1.26021 1.26021 

Voraussetzungen

keine Ausreisßer

Cook’s distance

  • Feste Grenzen
  • Werte < 1 sind unproblematisch

Leverage

  • Beziehen die Anzahl der Prädiktoren (k) und n mitein
  • (k + 1)/n
  • 2(k + 1)/n
  • 3(k + 1)/n

Voraussetzungen

keine Ausreisßer

Cook’s distance

Zeig mir den Code
plot(m1, 4, id.n = 5)

Voraussetzungen

… oder alles auf einmal

mit der check_model()-Funktion aus dem performance-Package

Zeig mir den Code
performance::check_model(m1)

Voraussetzungen

… oder alles auf einmal

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style G fill: #009e73

Ergebnis berichten

APA-konfrorm



Parameter

Fit

b

95% CI (b)

t

df

p

b*

95% CI (b*)

(Intercept)

12.42

[10.90, 13.94]

16.10

272

< .001***

0.00

[-0.10, 0.10]

hexp

-0.03

[-0.20, 0.15]

-0.30

272

.767

-0.02

[-0.13, 0.10]

relu

5.40

[4.28, 6.53]

9.48

272

< .001***

0.54

[0.43, 0.66]

AIC

1,631.27

AICc

1,631.42

BIC

1,645.74

R2

0.29

R2 (adj.)

0.28

Sigma

4.66

… jetzt mit Interaktion

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style G fill: #009e73

Multiple lineare Regression berechnen

Zentrierung der Prädiktoren

  • \(x_i - \bar{x}\)

Multiple lineare Regression berechnen

Zentrierung der Prädiktoren

Zeig mir den Code
# Exposition
ds$hexp_cen <- scale(ds$hexp,
                     scale = FALSE)

# Abneigung
ds$relu_cen <- scale(ds$relu,
                     scale = FALSE)

Multiple lineare Regression berechnen

Berechnen des Modells

Zeig mir den Code
m2 <- lm(mghs ~ hexp_cen * relu_cen, # "*" mit Interaktion
         data = ds)

Multiple lineare Regression berechnen

Berechnen des Modells

Was haben wir gefunden?

Zeig mir den Code
summary(m2)

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Graphisch

Zeig mir den Code
plot(m2, 1)

Voraussetzungen

Lineare Beziehung zwischen den Variablen

Statistisch - Rainbow-Test

Zeig mir den Code
lmtest::raintest(m2) 

    Rainbow test

data:  m2
Rain = 1.1316, df1 = 138, df2 = 133, p-value = 0.2368

Voraussetzungen

Normalverteilung der Residuen

Graphisch - Q-Q Plot

Zeig mir den Code
plot(m2, 2)

Voraussetzungen

Normalverteilung der Residuen

Statistisch - Shapiro-Wilk-Test

Zeig mir den Code
shapiro.test(m2$residuals)

    Shapiro-Wilk normality test

data:  m2$residuals
W = 0.99034, p-value = 0.06574

Voraussetzungen

Homoskedastizität

Graphisch

Zeig mir den Code
plot(m2, 3)

Voraussetzungen

Homoskedastizität

Statistisch - Breusch-Pagan-Test

Zeig mir den Code
car::ncvTest(m2)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.01463697, Df = 1, p = 0.9037

Voraussetzungen

Multikollinearität

Variance inflation factor (VIF)

Zeig mir den Code
car::vif(m2,
         type = "predictor")
         GVIF Df GVIF^(1/(2*Df)) Interacts With Other Predictors
hexp_cen    1  3               1       relu_cen             --  
relu_cen    1  3               1       hexp_cen             --  

Voraussetzungen

keine Ausreisßer

Cook’s distance

Zeig mir den Code
plot(m2, 4, id.n = 5)

Voraussestzungen

alles auf einmal

Interaktion

Simple slopes I

  • M
  • -1SD
  • 1SD
Zeig mir den Code
int <- interactions::probe_interaction(m2,
                                pred = hexp_cen, 
                                modx = relu_cen) 

int$simslopes

Interaktion

Simple slopes I

  • M
  • -1SD
  • 1SD
JOHNSON-NEYMAN INTERVAL 

When relu_cen is OUTSIDE the interval [-0.20, 0.08], the slope of hexp_cen
is p < .05.

Note: The range of observed values of relu_cen is [-1.22, 1.78]

SIMPLE SLOPES ANALYSIS 

Slope of hexp_cen when relu_cen = -5.536211e-01 (- 1 SD): 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -0.58   0.11    -5.47   0.00

Slope of hexp_cen when relu_cen = -2.014550e-16 (Mean): 

  Est.   S.E.   t val.      p
------ ------ -------- ------
  0.07   0.08     0.89   0.37

Slope of hexp_cen when relu_cen =  5.536211e-01 (+ 1 SD): 

  Est.   S.E.   t val.      p
------ ------ -------- ------
  0.72   0.12     5.89   0.00

Interaktion

Simple slopes I

  • M
  • -1SD
  • 1SD

Graphisch

Zeig mir den Code
int <- interactions::probe_interaction(m2,
                                pred = hexp_cen, 
                                modx = relu_cen) 

int$interactplot

Interaktion

Simple slopes I

Graphisch

Ablauf

flowchart TB

A(Daten einlesen) --> 
B(Übersicht verschaffen) --> 
C(Namen brauchbar machen) --> 
D(Ergeben die Daten Sinn?) --> 
E(Skalenwerte berechnen) --> 
F(multiple lineare Regression) --> 
G(Ergebnis berichten)

style G fill: #009e73

Ergebnis berichten

APA-konfrorm


Parameter

Fit

b

95% CI (b)

t

df

p

b*

95% CI (b*)

(Intercept)

12.42

[10.90, 13.94]

16.10

272

< .001***

0.00

[-0.10, 0.10]

hexp

-0.03

[-0.20, 0.15]

-0.30

272

.767

-0.02

[-0.13, 0.10]

relu

5.40

[4.28, 6.53]

9.48

272

< .001***

0.54

[0.43, 0.66]

AIC

1,631.27

AICc

1,631.42

BIC

1,645.74

R2

0.29

R2 (adj.)

0.28

Sigma

4.66